################################################################################################
# This script takes the bootstrapped QV data under the rules specified in the QV Welfare  
# Simulations Script and then produces graphs for some welfare measures as in 
# Figures 4 and 5 and Table 2 of the Columbia Elections Paper.
# 
# All graphs are exported to ./SV & QV Simulations (No Recalibration)/QV - Welfare & Minority Victories/
# 
# Author: Luis S.
# Last Modified: 7.25.16
################################################################################################

library(dplyr)
library(ggplot2)
library(ineq)
library(stargazer)

# setwd("C:/Users/Luis/Dropbox/Research/SV vs QV California Experiment/Experiment Analysis & Results/")

QV_EducBSPref <- read.csv("./SV & QV Simulations (No Recalibration)/QV - Raw Simulation Data/QV_EducBSPref.csv", header = T, stringsAsFactors = F)
QV_ImmBSPref <- read.csv("./SV & QV Simulations (No Recalibration)/QV - Raw Simulation Data/QV_ImmBSPref.csv", header = T, stringsAsFactors = F)
QV_BondBSPref <- read.csv("./SV & QV Simulations (No Recalibration)/QV - Raw Simulation Data/QV_BondBSPref.csv", header = T, stringsAsFactors = F)
QV_TeachBSPref <- read.csv("./SV & QV Simulations (No Recalibration)/QV - Raw Simulation Data/QV_TeachBSPref.csv", header = T, stringsAsFactors = F)

QV_PrefSummary_raw <- read.csv("./SV & QV Simulations (No Recalibration)/QV - Raw Simulation Data/QV_PrefPointSummary.csv", header = T, stringsAsFactors = F)

QV_RuleAResults <- read.csv("./SV & QV Simulations (No Recalibration)/QV - Raw Simulation Data/QV_RuleAResults.csv", header = T, stringsAsFactors = F)
QV_RuleBResults <- read.csv("./SV & QV Simulations (No Recalibration)/QV - Raw Simulation Data/QV_RuleBResults.csv", header = T, stringsAsFactors = F)
QV_RuleBResults_CNV <- read.csv("./SV & QV Simulations (No Recalibration)/QV - Raw Simulation Data/QV_RuleBResults_CNV.csv", header = T, stringsAsFactors = F)
QV_RuleCResults <- read.csv("./SV & QV Simulations (No Recalibration)/QV - Raw Simulation Data/QV_RuleCResults.csv", header = T, stringsAsFactors = F)
QV_RuleDResults_Unif <- read.csv("./SV & QV Simulations (No Recalibration)/QV - Raw Simulation Data/QV_RuleDResults_UniformVC.csv", header = T, stringsAsFactors = F)
QV_RuleDResults_Emp <- read.csv("./SV & QV Simulations (No Recalibration)/QV - Raw Simulation Data/QV_RuleDResults_EmpiricalVC.csv", header = T, stringsAsFactors = F)
QV_RuleDResults_5050 <- read.csv("SV & QV Simulations (No Recalibration)/QV - Raw Simulation Data/QV_RuleDResults_FiftyFifty.csv", header = T, stringsAsFactors = F)
QV_RuleEResults_Unif <- read.csv("./SV & QV Simulations (No Recalibration)/QV - Raw Simulation Data/QV_RuleEResults_UniformVC.csv", header = T, stringsAsFactors = F)
QV_RuleEResults_Emp <- read.csv("./SV & QV Simulations (No Recalibration)/QV - Raw Simulation Data/QV_RuleEResults_EmpiricalVC.csv", header = T, stringsAsFactors = F)

QV_RuleAResults <- mutate(QV_RuleAResults, Rule = "A")
QV_RuleBResults <- mutate(QV_RuleBResults, Rule = "B")
QV_RuleBResults_CNV <- mutate(QV_RuleBResults_CNV, Rule = "B - CNV")
QV_RuleCResults <- mutate(QV_RuleCResults, Rule = "C")
QV_RuleDResults_Unif <- mutate(QV_RuleDResults_Unif, Rule = "D - Unif")
QV_RuleDResults_Emp <- mutate(QV_RuleDResults_Emp, Rule = "D - Emp")
QV_RuleDResults_5050 <- mutate(QV_RuleDResults_5050, Rule = "D - 50/50")
QV_RuleEResults_Unif <- mutate(QV_RuleEResults_Unif, Rule = "E - Unif")
QV_RuleEResults_Emp <- mutate(QV_RuleEResults_Emp, Rule = "E - Emp")

QV_AllRules <- bind_rows(QV_RuleAResults, QV_RuleBResults, QV_RuleBResults_CNV, QV_RuleCResults,  QV_RuleDResults_Emp, QV_RuleDResults_Unif, QV_RuleDResults_5050, QV_RuleEResults_Emp, QV_RuleEResults_Unif)


# Print Average Votes and Points InFavor and Against across simulations for each issue
stargazer(as.data.frame(as.matrix(colMeans(QV_PrefSummary_raw[,-1]))) %>% tibble::rownames_to_column("Category") %>% rename(Mean = V1), 
          title = "Average of Votes and Points Per Issue (By Side)", type = "text", digits = 3, summary = FALSE, rownames = FALSE, 
          out = "./SV & QV Simulations (No Recalibration)/QV - Welfare & Minority Victories/Table - Average of Votes and Points Per Issue (By Side).txt")


####################################
###### Determine Minority Win ######
###################################%


QV_AllRules <- QV_AllRules %>% 
  mutate(QV1_MinorityWon = ifelse(QV1_Winner_BilingualEduc == "Minority" | QV1_Winner_Immigration == "Minority" | 
                                    QV1_Winner_PublicBonds == "Minority" | QV1_Winner_TeacherTenure == "Minority",T, F),
         QV2_MinorityWon = ifelse(QV2_Winner_BilingualEduc == "Minority" | QV2_Winner_Immigration == "Minority" | 
                                    QV2_Winner_PublicBonds == "Minority" | QV2_Winner_TeacherTenure == "Minority",T, F),
         MinorityWinIsEff = ifelse(NQV_Result_BilingualEduc != Eff_Result_BilingualEduc | NQV_Result_Immigration != Eff_Result_Immigration |
                                     NQV_Result_PublicBonds != Eff_Result_PublicBonds | NQV_Result_TeacherTenure != Eff_Result_TeacherTenure, T, F)) %>%
  mutate(QV1IsFullyEff = ifelse(QV1_Result_BilingualEduc != Eff_Result_BilingualEduc | QV1_Result_Immigration != Eff_Result_Immigration |
                                  QV1_Result_PublicBonds != Eff_Result_PublicBonds | QV1_Result_TeacherTenure != Eff_Result_TeacherTenure, F, T),
         QV2IsFullyEff = ifelse(QV2_Result_BilingualEduc != Eff_Result_BilingualEduc | QV2_Result_Immigration != Eff_Result_Immigration |
                                  QV2_Result_PublicBonds != Eff_Result_PublicBonds | QV2_Result_TeacherTenure != Eff_Result_TeacherTenure, F, T),
         QV1IsEff_BE = ifelse(QV1_Result_BilingualEduc == Eff_Result_BilingualEduc,T,F),
         QV1IsEff_IM = ifelse(QV1_Result_Immigration == Eff_Result_Immigration,T,F),
         QV1IsEff_PB = ifelse(QV1_Result_PublicBonds == Eff_Result_PublicBonds,T,F),
         QV1IsEff_TT = ifelse(QV1_Result_TeacherTenure == Eff_Result_TeacherTenure,T,F),
         QV2IsEff_BE = ifelse(QV2_Result_BilingualEduc == Eff_Result_BilingualEduc,T,F),
         QV2IsEff_IM = ifelse(QV2_Result_Immigration == Eff_Result_Immigration,T,F),
         QV2IsEff_PB = ifelse(QV2_Result_PublicBonds == Eff_Result_PublicBonds,T,F),
         QV2IsEff_TT = ifelse(QV2_Result_TeacherTenure == Eff_Result_TeacherTenure,T,F)) %>%
  mutate(QV1_MinorityWon_BE = ifelse(QV1_Winner_BilingualEduc == "Minority",T,F), 
         QV1_MinorityWon_IM = ifelse(QV1_Winner_Immigration == "Minority",T,F),
         QV1_MinorityWon_PB = ifelse(QV1_Winner_PublicBonds == "Minority",T,F),
         QV1_MinorityWon_TT = ifelse(QV1_Winner_TeacherTenure == "Minority",T, F),
         QV2_MinorityWon_BE = ifelse(QV2_Winner_BilingualEduc == "Minority",T,F), 
         QV2_MinorityWon_IM = ifelse(QV2_Winner_Immigration == "Minority",T,F),
         QV2_MinorityWon_PB = ifelse(QV2_Winner_PublicBonds == "Minority",T,F),
         QV2_MinorityWon_TT = ifelse(QV2_Winner_TeacherTenure == "Minority",T, F),
         MinorityWinIsEff_BE = ifelse(NQV_Result_BilingualEduc != Eff_Result_BilingualEduc,T,F), 
         MinorityWinIsEff_IM = ifelse(NQV_Result_Immigration != Eff_Result_Immigration,T,F), 
         MinorityWinIsEff_PB = ifelse(NQV_Result_PublicBonds != Eff_Result_PublicBonds,T,F), 
         MinorityWinIsEff_TT = ifelse(NQV_Result_TeacherTenure != Eff_Result_TeacherTenure, T, F)) %>%
  mutate(EffMinorityWinWithQV1_BE = ifelse(QV1_MinorityWon_BE & MinorityWinIsEff_BE, T, F),
         EffMinorityWinWithQV1_IM = ifelse(QV1_MinorityWon_IM & MinorityWinIsEff_IM, T, F),
         EffMinorityWinWithQV1_PB = ifelse(QV1_MinorityWon_PB & MinorityWinIsEff_PB, T, F),
         EffMinorityWinWithQV1_TT = ifelse(QV1_MinorityWon_TT & MinorityWinIsEff_TT, T, F),
         EffMinorityWinWithQV2_BE = ifelse(QV2_MinorityWon_BE & MinorityWinIsEff_BE, T, F),
         EffMinorityWinWithQV2_IM = ifelse(QV2_MinorityWon_IM & MinorityWinIsEff_IM, T, F),
         EffMinorityWinWithQV2_PB = ifelse(QV2_MinorityWon_PB & MinorityWinIsEff_PB, T, F),
         EffMinorityWinWithQV2_TT = ifelse(QV2_MinorityWon_TT & MinorityWinIsEff_TT, T, F)) %>%
  mutate(EffMinorityWinWithQV1 = ifelse(QV1_MinorityWon & MinorityWinIsEff, T, F), EffMinorityWinWithQV2 = ifelse(QV2_MinorityWon & MinorityWinIsEff, T, F))


#---- realized share of surplus (RSS)

# RSS_QV = (W_QV - R)/(W_Eff - R) where R = sum over all i,k preference intensities divided by 2 (to serve as a welfare floor)

welfFloor <- rowSums(QV_PrefSummary_raw %>% select(-matches("subjects")) %>% select(-matches("Round")), na.rm = TRUE)/2
welfFloor <- data.frame(Round = 1:10000, WelfFloor = welfFloor)

QV_AllRules <- full_join(QV_AllRules, welfFloor, by = "Round") 
QV_AllRules <- QV_AllRules %>% mutate(RSS_Maj = (AggregateWelfare_noQV - WelfFloor)/(AggregateWelfare_Eff - WelfFloor),
                                      RSS_QV1 = (AggregateWelfare_QV1 - WelfFloor)/(AggregateWelfare_Eff - WelfFloor),
                                      RSS_QV2 = (AggregateWelfare_QV2 - WelfFloor)/(AggregateWelfare_Eff - WelfFloor))




#---


QV_RuleAResults <- filter(QV_AllRules, Rule == "A")
QV_RuleBResults <- filter(QV_AllRules, Rule == "B")
QV_RuleBResults_CNV <- filter(QV_AllRules, Rule == "B - CNV")
QV_RuleCResults <- filter(QV_AllRules, Rule == "C")
QV_RuleDResults_Unif <- filter(QV_AllRules, Rule == "D - Unif")
QV_RuleDResults_Emp <- filter(QV_AllRules, Rule == "D - Emp")
QV_RuleDResults_5050 <- filter(QV_AllRules, Rule == "D - 50/50")
QV_RuleEResults_Unif <- filter(QV_AllRules, Rule == "E - Unif")
QV_RuleEResults_Emp <- filter(QV_AllRules, Rule == "E - Emp")

bootstrapIterations <- nrow(QV_RuleAResults)


######################################
##### Realized Share of Surplus ######
#####################################%

# QV

getNthEmpPercentile <- function(nums, NthPercent)
{
  orderedNums <- sort(nums)
  totalNums <- length(orderedNums)
  
  closestPercentile <- round((totalNums/100)*NthPercent)
  
  return(orderedNums[closestPercentile])
}


RSS_summary <- QV_AllRules %>% group_by(Rule) %>% 
  summarise(AvgRSS_Maj = mean(RSS_Maj), Maj_lower95CI = getNthEmpPercentile(RSS_Maj,2.5), Maj_upper95CI = getNthEmpPercentile(RSS_Maj,97.5),
            AvgRSS_QV1 = mean(RSS_QV1), QV1_lower95CI = getNthEmpPercentile(RSS_QV1,2.5), QV1_upper95CI = getNthEmpPercentile(RSS_QV1,97.5), 
            AvgRSS_QV2 = mean(RSS_QV2), QV2_lower95CI = getNthEmpPercentile(RSS_QV2,2.5), QV2_upper95CI = getNthEmpPercentile(RSS_QV2,97.5))


RSS_Test <- QV_AllRules %>% 
  mutate(Diff_QV1minusMaj = RSS_QV1 - RSS_Maj,
         Diff_QV2minusMaj = RSS_QV2 - RSS_Maj) %>%
  group_by(Rule) %>% 
  summarise(AvgDiff1_QV1minusMaj = mean(Diff_QV1minusMaj), Diff1_lower95CI = getNthEmpPercentile(Diff_QV1minusMaj,2.5), Diff1_upper95CI = getNthEmpPercentile(Diff_QV1minusMaj,97.5), 
            AvgDiff2_QV2minusMaj = mean(Diff_QV2minusMaj), Diff2_lower95CI = getNthEmpPercentile(Diff_QV2minusMaj,2.5), Diff2_upper95CI = getNthEmpPercentile(Diff_QV2minusMaj,97.5))




sink("./SV & QV Simulations (No Recalibration)/QV - Welfare & Minority Victories/Table - Realized Share of Surplus (By Rule).txt")
stargazer(as.data.frame(RSS_summary), title = "Realized Share of Surplus (RSS)", type = "text", digits = 4, summary = FALSE, rownames = FALSE)
stargazer(as.data.frame(RSS_Test), title = "Tests of Differences in RSS", type = "text", digits = 4, summary = FALSE, rownames = FALSE)
sink()


##################################################
####### Debug Table (per issue & per rule) #######
#################################################%

# 1. Number of simulations won by the InFavor side under majority
# 2. Number of sims won by the InFavor side under efficiency
# 3. Number of sims in which majority is efficient.
# 4. Number of sims won by the InFavor side under QV
# 5. Number of sims in which QV is efficient.
# 6. Number of sims in which QV is equivalent to majority.
# 7. Number of sims in which QV differs from maj and the InFavor side wins under QV.
# 8. Number of efficient minority victories under QV. 

ruleList <- c("A", "B", "B - CNV", "C", "D - Unif", "D - Emp", "D - 50/50","E - Unif", "E - Emp")
propList <- c("Immigration", "BilingualEduc", "PublicBonds", "TeacherTenure")
abbrevList <- c("IM", "BE", "PB", "TT")


finalResults <- data.frame()
# currentRule <- "A"

for(currentRule in ruleList)
{
  for(index in 1:4)
  {
    # index <- 1
    currentProp <- propList[index] 
    currentAbbrev <- abbrevList[index] 
    
    tempResults <- data.frame(Rule = currentRule, Proposal = currentProp, stringsAsFactors = F)
    
    # Number of simulations won by the InFavor side under majority
    temp <- QV_AllRules %>% filter(Rule == currentRule) %>% group_by_(paste("NQV_Result_",currentProp, sep = "")) %>%
      summarise(Count = n()) %>% filter_(paste("NQV_Result_",currentProp," == 'InFavor'", sep = "")) 
    try(tempResults$nSims_InF_Maj <- temp$Count, silent = T)
    
    # Number of sims won by the InFavor side under efficiency
    temp <- QV_AllRules %>% filter(Rule == currentRule) %>% group_by_(paste("Eff_Result_",currentProp, sep = "")) %>%
      summarise(Count = n()) %>% filter_(paste("Eff_Result_",currentProp," == 'InFavor'", sep = "")) 
    try(tempResults$nSims_InF_Eff <- temp$Count, silent = T)
    
    # Number of sims won by the InFavor side under, say, QV, rule A.
    temp <- QV_AllRules %>% filter(Rule == currentRule) %>% group_by_(paste("QV1_Result_",currentProp, sep = "")) %>%
      summarise(Count = n()) %>% filter_(paste("QV1_Result_",currentProp," == 'InFavor'", sep = "")) 
    try(tempResults$nSims_InF_QV1 <- temp$Count, silent = T)
    
    # Number of sims in which majority is efficient.
    temp <- QV_AllRules %>% filter(Rule == currentRule) %>% group_by_(paste("MinorityWinIsEff_",currentAbbrev, sep = "")) %>%
      summarise(Count = n()) %>% filter_(paste("MinorityWinIsEff_",currentAbbrev," == 'FALSE'", sep = "")) 
    try(tempResults$nSims_MajIsEff <- temp$Count, silent = T)
    
    # Number of sims in which QV is efficient.
    temp <- QV_AllRules %>% filter(Rule == currentRule) %>% group_by_(paste("QV1IsEff_",currentAbbrev, sep = "")) %>%
      summarise(Count = n()) %>% filter_(paste("QV1IsEff_",currentAbbrev," == 'TRUE'", sep = "")) 
    try(tempResults$nSims_QV1IsEff <- temp$Count, silent = T)
    
    # Number of sims in which QV is equivalent to majority.
    temp <- QV_AllRules %>% filter(Rule == currentRule) %>% 
      mutate(MajEqualsQV_BE = ifelse(NQV_Result_BilingualEduc == QV1_Result_BilingualEduc,T,F),
             MajEqualsQV_IM = ifelse(NQV_Result_Immigration == QV1_Result_Immigration,T,F),
             MajEqualsQV_PB = ifelse(NQV_Result_PublicBonds == QV1_Result_PublicBonds,T,F),
             MajEqualsQV_TT = ifelse(NQV_Result_TeacherTenure == QV1_Result_TeacherTenure,T,F)) %>%
      group_by_(paste("MajEqualsQV_",currentAbbrev, sep = "")) %>%
      summarise(Count = n()) %>% filter_(paste("MajEqualsQV_",currentAbbrev," == 'TRUE'", sep = "")) 
    try(tempResults$nSims_MajEqualsQV1 <- temp$Count, silent = T)
    
    # Number of sims in which QV differs from maj and the InFavor side wins under QV rule A.
    temp <- QV_AllRules %>% filter(Rule == currentRule) %>% 
      mutate(MajEqualsQV_BE = ifelse(NQV_Result_BilingualEduc == QV1_Result_BilingualEduc,T,F),
             MajEqualsQV_IM = ifelse(NQV_Result_Immigration == QV1_Result_Immigration,T,F),
             MajEqualsQV_PB = ifelse(NQV_Result_PublicBonds == QV1_Result_PublicBonds,T,F),
             MajEqualsQV_TT = ifelse(NQV_Result_TeacherTenure == QV1_Result_TeacherTenure,T,F)) %>%
      group_by_(paste("MajEqualsQV_",currentAbbrev, sep = ""), paste("QV1_Result_",currentProp, sep = "")) %>%
      summarise(Count = n()) %>% filter_(paste("MajEqualsQV_",currentAbbrev," == 'FALSE'", sep = ""), 
                                         paste("QV1_Result_",currentProp," == 'InFavor'", sep = "")) 
    try(tempResults$nSims_MajDiffQV1_QVInF <- temp$Count, silent = T)
    
    
    # Number of efficient minority victories under QV rule A. 
    temp <- QV_AllRules %>% filter(Rule == currentRule) %>% group_by_(paste("EffMinorityWinWithQV1_",currentAbbrev, sep = "")) %>%
      summarise(Count = n())%>% filter_(paste("EffMinorityWinWithQV1_",currentAbbrev," == 'TRUE'", sep = "")) 
    try(tempResults$nSims_EffMinWinQV1 <- temp$Count, silent = T)
    
    # Bind results
    finalResults <- bind_rows(finalResults, tempResults)
  }
}

finalResults[is.na(finalResults)] <- 0

write.csv(finalResults, "SV & QV Simulations (No Recalibration)/QV - Welfare & Minority Victories/Table - Summary of Sims (per issue & per rule).csv", row.names = F)


##############################################################
####### Prop of Victories In Favor of Issue (per rule) #######
#############################################################%

# NQV

temp1 <- QV_AllRules %>% group_by(Rule, NQV_Result_Immigration) %>% summarise(Count = n()) %>% 
  mutate(Prop = Count/sum(Count)) %>% filter(NQV_Result_Immigration == "InFavor") %>%
  select(Rule, Prop) %>% rename(IM = Prop)

temp2 <- QV_AllRules %>% group_by(Rule, NQV_Result_BilingualEduc) %>% summarise(Count = n()) %>% 
  mutate(Prop = Count/sum(Count)) %>% filter(NQV_Result_BilingualEduc == "InFavor") %>%
  select(Rule, Prop) %>% rename(BE = Prop)

temp3 <- QV_AllRules %>% group_by(Rule, NQV_Result_PublicBonds) %>% summarise(Count = n()) %>% 
  mutate(Prop = Count/sum(Count)) %>% filter(NQV_Result_PublicBonds == "InFavor") %>%
  select(Rule, Prop) %>% rename(PB = Prop)

temp4 <- QV_AllRules %>% group_by(Rule, NQV_Result_TeacherTenure) %>% summarise(Count = n()) %>% 
  mutate(Prop = Count/sum(Count)) %>% filter(NQV_Result_TeacherTenure == "InFavor") %>%
  select(Rule, Prop) %>% rename(TT = Prop)

SimpleMajSummary <- full_join(full_join(temp1,temp2, by = "Rule"),full_join(temp3,temp4, by = "Rule"), by = "Rule")
SimpleMajSummary <- SimpleMajSummary %>% mutate(Vote = "Simple Majority")
SimpleMajSummary[is.na(SimpleMajSummary)] <- 0

# QV1

temp1 <- QV_AllRules %>% group_by(Rule, QV1_Result_Immigration) %>% summarise(Count = n()) %>% 
  mutate(Prop = Count/sum(Count)) %>% filter(QV1_Result_Immigration == "InFavor") %>%
  select(Rule, Prop) %>% rename(IM = Prop)

temp2 <- QV_AllRules %>% group_by(Rule, QV1_Result_BilingualEduc) %>% summarise(Count = n()) %>% 
  mutate(Prop = Count/sum(Count)) %>% filter(QV1_Result_BilingualEduc == "InFavor") %>%
  select(Rule, Prop) %>% rename(BE = Prop)

temp3 <- QV_AllRules %>% group_by(Rule, QV1_Result_PublicBonds) %>% summarise(Count = n()) %>% 
  mutate(Prop = Count/sum(Count)) %>% filter(QV1_Result_PublicBonds == "InFavor") %>%
  select(Rule, Prop) %>% rename(PB = Prop)

temp4 <- QV_AllRules %>% group_by(Rule, QV1_Result_TeacherTenure) %>% summarise(Count = n()) %>% 
  mutate(Prop = Count/sum(Count)) %>% filter(QV1_Result_TeacherTenure == "InFavor") %>%
  select(Rule, Prop) %>% rename(TT = Prop)

Vote1Summary <- full_join(full_join(temp1,temp2, by = "Rule"),full_join(temp3,temp4, by = "Rule"), by = "Rule")
Vote1Summary <- Vote1Summary %>% mutate(Vote = "QV1")
Vote1Summary[is.na(Vote1Summary)] <- 0


# QV2

temp1 <- QV_AllRules %>% group_by(Rule, QV2_Result_Immigration) %>% summarise(Count = n()) %>% 
  mutate(Prop = Count/sum(Count)) %>% filter(QV2_Result_Immigration == "InFavor") %>%
  select(Rule, Prop) %>% rename(IM = Prop)

temp2 <- QV_AllRules %>% group_by(Rule, QV2_Result_BilingualEduc) %>% summarise(Count = n()) %>% 
  mutate(Prop = Count/sum(Count)) %>% filter(QV2_Result_BilingualEduc == "InFavor") %>%
  select(Rule, Prop) %>% rename(BE = Prop)

temp3 <- QV_AllRules %>% group_by(Rule, QV2_Result_PublicBonds) %>% summarise(Count = n()) %>% 
  mutate(Prop = Count/sum(Count)) %>% filter(QV2_Result_PublicBonds == "InFavor") %>%
  select(Rule, Prop) %>% rename(PB = Prop)

temp4 <- QV_AllRules %>% group_by(Rule, QV2_Result_TeacherTenure) %>% summarise(Count = n()) %>% 
  mutate(Prop = Count/sum(Count)) %>% filter(QV2_Result_TeacherTenure == "InFavor") %>%
  select(Rule, Prop) %>% rename(TT = Prop)

Vote2Summary <- full_join(full_join(temp1,temp2, by = "Rule"),full_join(temp3,temp4, by = "Rule"), by = "Rule")
Vote2Summary <- Vote2Summary %>% mutate(Vote = "QV2")
Vote2Summary[is.na(Vote2Summary)] <- 0

# Combine and Print

sink("SV & QV Simulations (No Recalibration)/QV - Welfare & Minority Victories/Table - Prop of Victories In Favor of Each Proposal (per Rule & Vote).txt")
stargazer(as.data.frame(SimpleMajSummary), type = "text", summary = FALSE, rownames = FALSE, digits = 3)
stargazer(as.data.frame(Vote1Summary), type = "text", summary = FALSE, rownames = FALSE, digits = 3)
stargazer(as.data.frame(Vote2Summary), type = "text", summary = FALSE, rownames = FALSE, digits = 3)
sink()


##############################################
####### Freq Histogram of Minority Win #######
#############################################%

#### Frequency of elections in which at least one minority victory occurs per rule (Fig 4 in Elections paper). ####%

# QV1
QV1_AllRules <- group_by(QV_AllRules, Rule, QV1_MinorityWon) %>%
  summarize(Freq = round(n()/bootstrapIterations, digits = 3)*100) %>%
  filter(QV1_MinorityWon == TRUE)

welfarePlot <- ggplot(data = QV1_AllRules, aes(x=Rule, y = Freq)) + 
  geom_bar(stat="identity") + scale_y_continuous(limits = c(0,100), breaks = seq(0, 100, by = 10)) +  
  xlab("Bonus Vote - Rules") + ylab("Frequency of Minority Victories") + 
  theme(axis.text.x=element_text(size=12, vjust = -.05), legend.text=element_text(size=16), title=element_text(size=12)) +
  ggtitle("QV1 (No Recal) - Bootstrapped Freq. of Minority Victories")
print(welfarePlot)
ggsave(welfarePlot, file="./SV & QV Simulations (No Recalibration)/QV - Welfare & Minority Victories/QV1 (No Recal) - Freq. Minority Victories.png", height = 5.5, width = 7.5)


# QV2
QV2_AllRules <- group_by(QV_AllRules, Rule, QV2_MinorityWon) %>%
  summarize(Freq = round(n()/bootstrapIterations, digits = 3)*100) %>%
  filter(QV2_MinorityWon == TRUE)

welfarePlot <- ggplot(data = QV2_AllRules, aes(x=Rule, y = Freq)) + 
  geom_bar(stat="identity") + scale_y_continuous(limits = c(0,100), breaks = seq(0, 100, by = 10)) +  
  xlab("Bonus Vote - Rules") + ylab("Frequency of Minority Victories") + 
  theme(axis.text.x=element_text(size=12, vjust = -.05), legend.text=element_text(size=16), title=element_text(size=12)) +
  ggtitle("QV2 (No Recal) - Bootstrapped Freq. of Minority Victories")
print(welfarePlot)
ggsave(welfarePlot, file="./SV & QV Simulations (No Recalibration)/QV - Welfare & Minority Victories/QV2 (No Recal) - Freq. Minority Victories.png", height = 5.5, width = 7.5)


###################################################
####### Freq Table of Minority Win (byRule) #######
##################################################%

#### Frequency Table of elections in which at least one minority victory occurs per rule (Table 2 in Elections paper). ###%

# QV1
Summary_Immigration <- group_by(QV_AllRules, Rule, QV1_Winner_Immigration) %>%
  summarize(QV1_Immigration_Freq = round(n()/bootstrapIterations, digits = 4)*100) %>% 
  filter(QV1_Winner_Immigration == "Majority") %>% mutate(QV1_Immigration_Freq = 100 - QV1_Immigration_Freq)

Summary_PublicBonds <- group_by(QV_AllRules, Rule, QV1_Winner_PublicBonds) %>%
  summarize(QV1_PublicBonds_Freq = round(n()/bootstrapIterations, digits = 4)*100) %>% 
  filter(QV1_Winner_PublicBonds == "Majority") %>% mutate(QV1_PublicBonds_Freq = 100 - QV1_PublicBonds_Freq)

Summary_TeacherTenure <- group_by(QV_AllRules, Rule, QV1_Winner_TeacherTenure) %>%
  summarize(QV1_TeacherTenure_Freq = round(n()/bootstrapIterations, digits = 4)*100) %>% 
  filter(QV1_Winner_TeacherTenure == "Majority") %>% mutate(QV1_TeacherTenure_Freq = 100 - QV1_TeacherTenure_Freq)

Summary_BilingualEduc <- group_by(QV_AllRules, Rule, QV1_Winner_BilingualEduc) %>%
  summarize(QV1_BilingualEduc_Freq = round(n()/bootstrapIterations, digits = 4)*100) %>% 
  filter(QV1_Winner_BilingualEduc == "Majority") %>% mutate(QV1_BilingualEduc_Freq = 100 - QV1_BilingualEduc_Freq)

temp1 <- full_join(Summary_Immigration, Summary_BilingualEduc, by = "Rule")
temp2 <- full_join(Summary_TeacherTenure, Summary_PublicBonds, by = "Rule")

QV1_MinorityVictoryRules <- full_join(temp1, temp2, by = "Rule") %>% 
  select(Rule, matches("Freq"))

stargazer(as.data.frame(QV1_MinorityVictoryRules), type = "text", digits = 3, summary = FALSE, rownames = FALSE)

# QV2
Summary_Immigration <- group_by(QV_AllRules, Rule, QV2_Winner_Immigration) %>%
  summarize(QV2_Immigration_Freq = round(n()/bootstrapIterations, digits = 4)*100) %>% 
  filter(QV2_Winner_Immigration == "Majority") %>% mutate(QV2_Immigration_Freq = 100 - QV2_Immigration_Freq)

Summary_PublicBonds <- group_by(QV_AllRules, Rule, QV2_Winner_PublicBonds) %>%
  summarize(QV2_PublicBonds_Freq = round(n()/bootstrapIterations, digits = 4)*100) %>% 
  filter(QV2_Winner_PublicBonds == "Majority") %>% mutate(QV2_PublicBonds_Freq = 100 - QV2_PublicBonds_Freq)

Summary_TeacherTenure <- group_by(QV_AllRules, Rule, QV2_Winner_TeacherTenure) %>%
  summarize(QV2_TeacherTenure_Freq = round(n()/bootstrapIterations, digits = 4)*100) %>% 
  filter(QV2_Winner_TeacherTenure == "Majority") %>% mutate(QV2_TeacherTenure_Freq = 100 - QV2_TeacherTenure_Freq)

Summary_BilingualEduc <- group_by(QV_AllRules, Rule, QV2_Winner_BilingualEduc) %>%
  summarize(QV2_BilingualEduc_Freq = round(n()/bootstrapIterations, digits = 4)*100) %>% 
  filter(QV2_Winner_BilingualEduc == "Majority") %>% mutate(QV2_BilingualEduc_Freq = 100 - QV2_BilingualEduc_Freq)

temp1 <- full_join(Summary_Immigration, Summary_BilingualEduc, by = "Rule")
temp2 <- full_join(Summary_TeacherTenure, Summary_PublicBonds, by = "Rule")

QV2_MinorityVictoryRules <- full_join(temp1, temp2, by = "Rule") %>% 
  select(Rule, matches("Freq"))

stargazer(as.data.frame(QV2_MinorityVictoryRules), type = "text", digits = 3, summary = FALSE, rownames = FALSE)

# Print them out to the folder 
sink("./SV & QV Simulations (No Recalibration)/QV - Welfare & Minority Victories/Table - Percent of Minority Victories (per Issue).txt")
stargazer(as.data.frame(QV1_MinorityVictoryRules), type = "text", summary = FALSE, rownames = FALSE, digits = 3)
stargazer(as.data.frame(QV2_MinorityVictoryRules), type = "text", summary = FALSE, rownames = FALSE, digits = 3)
sink()







######################################################################
####### Freq Table of Efficient Minority Wins (byRule & Issue) #######
#####################################################################%

# We want two things:
# - proportion of times in which a minority win would have been efficient
# - proportion of times in which a minority win would have been efficient AND was achieved by BV

#--- Prop of Times a Minority Win is Efficient
temp_BE <- QV_RuleAResults %>% group_by(MinorityWinIsEff_BE) %>% summarise(BE_Count = n()) %>% rename(MinorityWinIsEff = MinorityWinIsEff_BE)
temp_IM <- QV_RuleAResults %>% group_by(MinorityWinIsEff_IM) %>% summarise(IM_Count = n()) %>% rename(MinorityWinIsEff = MinorityWinIsEff_IM)
temp_PB <- QV_RuleAResults %>% group_by(MinorityWinIsEff_PB) %>% summarise(PB_Count = n()) %>% rename(MinorityWinIsEff = MinorityWinIsEff_PB)
temp_TT <- QV_RuleAResults %>% group_by(MinorityWinIsEff_TT) %>% summarise(TT_Count = n()) %>% rename(MinorityWinIsEff = MinorityWinIsEff_TT)

MinorityWinEff <- QV_RuleAResults %>% group_by(MinorityWinIsEff) %>% summarise(Count = n()) %>% mutate(Prop = Count/sum(Count)) %>%
  full_join(temp_BE, by = "MinorityWinIsEff") %>% full_join(temp_IM, by = "MinorityWinIsEff") %>%
  full_join(temp_PB, by = "MinorityWinIsEff") %>% full_join(temp_TT, by = "MinorityWinIsEff")
MinorityWinEff[is.na(MinorityWinEff)] <- 0

#--- Prop of Times a Minority Win is achieved by QV

# QV1
temp_BE <- QV_AllRules %>% group_by(Rule, QV1_MinorityWon_BE) %>% summarise(BE_Count = n()) %>% rename(MinorityWinWithQV = QV1_MinorityWon_BE)
temp_IM <- QV_AllRules %>% group_by(Rule, QV1_MinorityWon_IM) %>% summarise(IM_Count = n()) %>% rename(MinorityWinWithQV = QV1_MinorityWon_IM)
temp_PB <- QV_AllRules %>% group_by(Rule, QV1_MinorityWon_PB) %>% summarise(PB_Count = n()) %>% rename(MinorityWinWithQV = QV1_MinorityWon_PB)
temp_TT <- QV_AllRules %>% group_by(Rule, QV1_MinorityWon_TT) %>% summarise(TT_Count = n()) %>% rename(MinorityWinWithQV = QV1_MinorityWon_TT)

MinorityWinByQV1 <- QV_AllRules %>% group_by(Rule, QV1_MinorityWon) %>% summarise(Count = n()) %>%
  mutate(Prop = Count/sum(Count), Vote = "QV1") %>% rename(MinorityWinWithQV = QV1_MinorityWon) %>%
  full_join(temp_BE, by = c("MinorityWinWithQV" = "MinorityWinWithQV", "Rule" = "Rule")) %>% 
  full_join(temp_IM, by = c("MinorityWinWithQV" = "MinorityWinWithQV", "Rule" = "Rule")) %>%
  full_join(temp_PB, by = c("MinorityWinWithQV" = "MinorityWinWithQV", "Rule" = "Rule")) %>% 
  full_join(temp_TT, by = c("MinorityWinWithQV" = "MinorityWinWithQV", "Rule" = "Rule"))


# QV2

temp_BE <- QV_AllRules %>% group_by(Rule, QV2_MinorityWon_BE) %>% summarise(BE_Count = n()) %>% rename(MinorityWinWithQV = QV2_MinorityWon_BE)
temp_IM <- QV_AllRules %>% group_by(Rule, QV2_MinorityWon_IM) %>% summarise(IM_Count = n()) %>% rename(MinorityWinWithQV = QV2_MinorityWon_IM)
temp_PB <- QV_AllRules %>% group_by(Rule, QV2_MinorityWon_PB) %>% summarise(PB_Count = n()) %>% rename(MinorityWinWithQV = QV2_MinorityWon_PB)
temp_TT <- QV_AllRules %>% group_by(Rule, QV2_MinorityWon_TT) %>% summarise(TT_Count = n()) %>% rename(MinorityWinWithQV = QV2_MinorityWon_TT)

MinorityWinByQV2 <- QV_AllRules %>% group_by(Rule, QV2_MinorityWon) %>% summarise(Count = n()) %>%
  mutate(Prop = Count/sum(Count), Vote = "QV2") %>% rename(MinorityWinWithQV = QV2_MinorityWon) %>%
  full_join(temp_BE, by = c("MinorityWinWithQV" = "MinorityWinWithQV", "Rule" = "Rule")) %>% 
  full_join(temp_IM, by = c("MinorityWinWithQV" = "MinorityWinWithQV", "Rule" = "Rule")) %>%
  full_join(temp_PB, by = c("MinorityWinWithQV" = "MinorityWinWithQV", "Rule" = "Rule")) %>% 
  full_join(temp_TT, by = c("MinorityWinWithQV" = "MinorityWinWithQV", "Rule" = "Rule"))

# Merged 

MinorityWinByQV <- bind_rows(MinorityWinByQV1, MinorityWinByQV2) %>% filter(MinorityWinWithQV == T)
MinorityWinByQV[is.na(MinorityWinByQV)] <- 0

#--- Prop of Times a Minority Win is Efficient AND was achieved by QV

# QV1

temp_BE <- QV_AllRules %>% group_by(Rule, EffMinorityWinWithQV1_BE) %>% summarise(BE_Count = n()) %>% rename(EffMinorityWinWithQV = EffMinorityWinWithQV1_BE)
temp_IM <- QV_AllRules %>% group_by(Rule, EffMinorityWinWithQV1_IM) %>% summarise(IM_Count = n()) %>% rename(EffMinorityWinWithQV = EffMinorityWinWithQV1_IM)
temp_PB <- QV_AllRules %>% group_by(Rule, EffMinorityWinWithQV1_PB) %>% summarise(PB_Count = n()) %>% rename(EffMinorityWinWithQV = EffMinorityWinWithQV1_PB)
temp_TT <- QV_AllRules %>% group_by(Rule, EffMinorityWinWithQV1_TT) %>% summarise(TT_Count = n()) %>% rename(EffMinorityWinWithQV = EffMinorityWinWithQV1_TT)

EffMinorityWinByQV1 <- QV_AllRules %>% group_by(Rule, EffMinorityWinWithQV1) %>% summarise(Count = n()) %>%
  mutate(Prop = Count/sum(Count), Vote = "QV1") %>% rename(EffMinorityWinWithQV = EffMinorityWinWithQV1)  %>%
  full_join(temp_BE, by = c("EffMinorityWinWithQV" = "EffMinorityWinWithQV", "Rule" = "Rule")) %>% 
  full_join(temp_IM, by = c("EffMinorityWinWithQV" = "EffMinorityWinWithQV", "Rule" = "Rule")) %>%
  full_join(temp_PB, by = c("EffMinorityWinWithQV" = "EffMinorityWinWithQV", "Rule" = "Rule")) %>% 
  full_join(temp_TT, by = c("EffMinorityWinWithQV" = "EffMinorityWinWithQV", "Rule" = "Rule"))

# QV2

temp_BE <- QV_AllRules %>% group_by(Rule, EffMinorityWinWithQV2_BE) %>% summarise(BE_Count = n()) %>% rename(EffMinorityWinWithQV = EffMinorityWinWithQV2_BE)
temp_IM <- QV_AllRules %>% group_by(Rule, EffMinorityWinWithQV2_IM) %>% summarise(IM_Count = n()) %>% rename(EffMinorityWinWithQV = EffMinorityWinWithQV2_IM)
temp_PB <- QV_AllRules %>% group_by(Rule, EffMinorityWinWithQV2_PB) %>% summarise(PB_Count = n()) %>% rename(EffMinorityWinWithQV = EffMinorityWinWithQV2_PB)
temp_TT <- QV_AllRules %>% group_by(Rule, EffMinorityWinWithQV2_TT) %>% summarise(TT_Count = n()) %>% rename(EffMinorityWinWithQV = EffMinorityWinWithQV2_TT)

EffMinorityWinByQV2 <- QV_AllRules %>% group_by(Rule, EffMinorityWinWithQV2) %>% summarise(Count = n()) %>%
  mutate(Prop = Count/sum(Count), Vote = "QV2") %>% rename(EffMinorityWinWithQV = EffMinorityWinWithQV2)  %>%
  full_join(temp_BE, by = c("EffMinorityWinWithQV" = "EffMinorityWinWithQV", "Rule" = "Rule")) %>% 
  full_join(temp_IM, by = c("EffMinorityWinWithQV" = "EffMinorityWinWithQV", "Rule" = "Rule")) %>%
  full_join(temp_PB, by = c("EffMinorityWinWithQV" = "EffMinorityWinWithQV", "Rule" = "Rule")) %>% 
  full_join(temp_TT, by = c("EffMinorityWinWithQV" = "EffMinorityWinWithQV", "Rule" = "Rule"))

# Merged 

EffMinorityWinByQV <- bind_rows(EffMinorityWinByQV1, EffMinorityWinByQV2) %>% filter(EffMinorityWinWithQV == T)
EffMinorityWinByQV[is.na(EffMinorityWinByQV)] <- 0

# Print them out to the folder
sink("SV & QV Simulations (No Recalibration)/QV - Welfare & Minority Victories/Table - Freq of Efficient Minority Victories.txt")
stargazer(as.data.frame(MinorityWinEff), title = "Prop of times in which Min win is efficient", type = "text", summary = FALSE, rownames = FALSE, digits = 3)
stargazer(as.data.frame(MinorityWinByQV), title = "Prop of times in which QV leads the minority to win (need not be eff)", type = "text", summary = FALSE, rownames = FALSE, digits = 3)
stargazer(as.data.frame(EffMinorityWinByQV), title = "Prop of times in which Min win is efficient and was achieved with QV", type = "text", summary = FALSE, rownames = FALSE, digits = 3)
sink()






######################################
####### Potential Welfare Gain #######
#####################################%

bins25 <- data.frame(Order = seq(0,25,1))
bins25 <- bins25 %>% mutate(Bins = cut(Order, seq(0,25,1), include.lowest = F))

# Potential Welfare Gain: it is the same for all rules, only rule A is used
# Defined as [W(eff) - w(maj)]/w(maj) for all sims
PotentialWelfareGain <- QV_RuleAResults %>%  
  mutate(WelfareFrontier = (AggregateWelfare_Eff - AggregateWelfare_noQV)/AggregateWelfare_noQV) %>%
  select(Round, AggregateWelfare_Eff, AggregateWelfare_noQV, WelfareFrontier) %>%
  mutate(WelfareFrontier = round(WelfareFrontier*100, digits = 2)) %>%
  mutate(Bins = cut(WelfareFrontier, seq(0,25,1), include.lowest = F)) %>%
  group_by(Bins) %>% summarise(Count = n()) %>% mutate(Freq = Count/sum(Count)) 

UnbinnedPotentialWelfareGain <- QV_RuleAResults %>%  
  mutate(WelfareFrontier = (AggregateWelfare_Eff - AggregateWelfare_noQV)/AggregateWelfare_noQV) %>%
  select(Round, AggregateWelfare_Eff, AggregateWelfare_noQV, WelfareFrontier) %>%
  mutate(WelfareFrontier = round(WelfareFrontier*100, digits = 2)) %>%
  mutate(Bins = cut(WelfareFrontier, seq(0,25,1), include.lowest = F))

l1 <- round(mean(na.omit(UnbinnedPotentialWelfareGain)$WelfareFrontier), digits = 2)
l2 <- round(median(na.omit(UnbinnedPotentialWelfareGain)$WelfareFrontier), digits = 2)
l1 <- paste("Mean PWG (% Maj. Welfare): ", l1, sep = "")
l2 <- paste("Median PWG (% Maj. Welfare): ", l2, sep = "")

temp <- left_join(bins25, PotentialWelfareGain, by = "Bins") %>% filter(Order != 0)
temp[is.na(temp)] <- 0

welfarePlot <- ggplot(data = na.omit(PotentialWelfareGain), aes(x=Bins, y = Freq)) + geom_bar(stat = "identity") +
  ggtitle("Relative Freq. of Potential Welfare Gain as % of Simple Majority Welfare (No Recal)") + 
  xlab("Potential Welfare Gain (as % of maj. welfare)") + ylab("Relative Frequency") +
  theme(axis.text.x=element_text(size=12, vjust = -.05, angle = 270), title=element_text(size=12)) + 
  annotate("text", x = 18, y = c(.1,.09), label = c(l1,l2)) + coord_cartesian(ylim = c(0, .12)) + 
  scale_x_discrete(drop=FALSE)
print(welfarePlot)
ggsave(welfarePlot, file="./SV & QV Simulations (No Recalibration)/QV - Welfare & Minority Victories/Potential Welfare Gain as Fraction of Simple Majority Welfare (No Recal).png", height = 5.5, width = 9)


######################################
####### Realized Welfare Gains #######
#####################################%


# (a) Majority: W(M)/W(eff) for each simulation. And average over all simulations.
# (b) QV: W(QV)/W(eff) for each simulation. Then both average over all simulations, and average over all s(m, QV) only.

#--- Modified ---#

getNthEmpPercentile <- function(nums, NthPercent)
{
  orderedNums <- sort(nums)
  totalNums <- length(orderedNums)
  
  closestPercentile <- round((totalNums/100)*NthPercent)
  
  return(orderedNums[closestPercentile])
}

NQV_RealizedWelfareGain <- QV_RuleAResults %>%
  mutate(RealWelfGain_All = AggregateWelfare_noQV/AggregateWelfare_Eff) %>%
  summarise(MeanRealWelfGain_All = mean(RealWelfGain_All), 
            All_lower95CI = getNthEmpPercentile(RealWelfGain_All,2.5),
            All_upper95CI = getNthEmpPercentile(RealWelfGain_All,97.5)) %>%
  mutate(Vote = "Simple Maj.")


QV1_RealizedWelfareGain <- QV_AllRules %>%
  mutate(RealWelfGain_All = AggregateWelfare_QV1/AggregateWelfare_Eff,
         RealWelfGain_MinQV = ifelse(QV1_MinorityWon == T, AggregateWelfare_QV1/AggregateWelfare_Eff, NA))  %>%
  group_by(Rule) %>% 
  summarise(MeanRealWelfGain_All = mean(RealWelfGain_All), 
            All_lower95CI = getNthEmpPercentile(RealWelfGain_All,2.5),
            All_upper95CI = getNthEmpPercentile(RealWelfGain_All,97.5),
            MeanRealWelfGain_MinQV = mean(RealWelfGain_MinQV, na.rm = T), 
            MinQV_lower95CI = getNthEmpPercentile(RealWelfGain_MinQV,2.5),
            MinQV_upper95CI = getNthEmpPercentile(RealWelfGain_MinQV,97.5)) %>%
  mutate(Vote = "QV1")


QV2_RealizedWelfareGain <- QV_AllRules %>%
  mutate(RealWelfGain_All = AggregateWelfare_QV2/AggregateWelfare_Eff,
         RealWelfGain_MinQV = ifelse(QV2_MinorityWon == T, AggregateWelfare_QV2/AggregateWelfare_Eff, NA))  %>%
  group_by(Rule) %>% 
  summarise(MeanRealWelfGain_All = mean(RealWelfGain_All), 
            All_lower95CI = getNthEmpPercentile(RealWelfGain_All,2.5),
            All_upper95CI = getNthEmpPercentile(RealWelfGain_All,97.5),
            MeanRealWelfGain_MinQV = mean(RealWelfGain_MinQV, na.rm = T), 
            MinQV_lower95CI = getNthEmpPercentile(RealWelfGain_MinQV,2.5),
            MinQV_upper95CI = getNthEmpPercentile(RealWelfGain_MinQV,97.5)) %>%
  mutate(Vote = "QV2")


AvgRealizedWelfareGain <- bind_rows(QV1_RealizedWelfareGain, QV2_RealizedWelfareGain, NQV_RealizedWelfareGain)


sink("./SV & QV Simulations (No Recalibration)/QV - Welfare & Minority Victories/Table - Average Realized Welfare Gain (per Rule).txt")
stargazer(as.data.frame(AvgRealizedWelfareGain), type = "text", summary = FALSE, rownames = FALSE, digits = 5)
sink()

AvgForegoneWelfareGain <- AvgRealizedWelfareGain %>% 
  mutate(MeanForegoneGain_All = 1 - MeanRealWelfGain_All, MeanForegoneGain_MinQV = 1 - MeanRealWelfGain_MinQV) %>%
  select(-matches("Real")) %>% select(-matches("CI"))


sink("./SV & QV Simulations (No Recalibration)/QV - Welfare & Minority Victories/Table - Average Foregone Gain (per Rule).txt")
stargazer(as.data.frame(AvgForegoneWelfareGain), type = "text", summary = FALSE, rownames = FALSE, digits = 5)
sink()


#--- Modified (end) ---#

# NQV_RealizedWelfareGain <- QV_RuleAResults %>%
#   mutate(RealWelfGain_All = AggregateWelfare_noQV/AggregateWelfare_Eff) %>%
#   summarise(MeanRealWelfGain_All = mean(RealWelfGain_All)) %>%
#   mutate(Vote = "Simple Maj.")
# 
# QV1_RealizedWelfareGain <- QV_AllRules %>%
#   mutate(RealWelfGain_All = AggregateWelfare_QV1/AggregateWelfare_Eff,
#          RealWelfGain_MinQV = ifelse(QV1_MinorityWon == T, AggregateWelfare_QV1/AggregateWelfare_Eff, NA))  %>%
#   group_by(Rule) %>% 
#   summarise(MeanRealWelfGain_All = mean(RealWelfGain_All),
#             MeanRealWelfGain_MinQV = mean(RealWelfGain_MinQV, na.rm = T)) %>%
#   mutate(Vote = "QV1")
# 
# QV2_RealizedWelfareGain <- QV_AllRules %>%
#   mutate(RealWelfGain_All = AggregateWelfare_QV2/AggregateWelfare_Eff,
#          RealWelfGain_MinQV = ifelse(QV2_MinorityWon == T, AggregateWelfare_QV2/AggregateWelfare_Eff, NA))  %>%
#   group_by(Rule) %>% 
#   summarise(MeanRealWelfGain_All = mean(RealWelfGain_All),
#             MeanRealWelfGain_MinQV = mean(RealWelfGain_MinQV, na.rm = T)) %>%
#   mutate(Vote = "QV2")
# 
# 
# AvgRealizedWelfareGain <- bind_rows(QV1_RealizedWelfareGain, QV2_RealizedWelfareGain, NQV_RealizedWelfareGain)
# 
# 
# 
# sink("./SV & QV Simulations (No Recalibration)/QV - Welfare & Minority Victories/Table - Average Realized Welfare Gain (per Rule).txt")
# stargazer(as.data.frame(AvgRealizedWelfareGain), type = "text", summary = FALSE, rownames = FALSE, digits = 5)
# sink()
# 
# AvgForegoneWelfareGain <- AvgRealizedWelfareGain %>% 
#   mutate(MeanForegoneGain_All = 1 - MeanRealWelfGain_All, MeanForegoneGain_MinSV = 1 - MeanRealWelfGain_MinQV) %>%
#   select(-matches("Real"))
# 
# 
# sink("./SV & QV Simulations (No Recalibration)/QV - Welfare & Minority Victories/Table - Average Foregone Gain (per Rule).txt")
# stargazer(as.data.frame(AvgForegoneWelfareGain), type = "text", summary = FALSE, rownames = FALSE, digits = 5)
# sink()




###################################################
####### Freq. of Welfare Increases (byRule) #######
##################################################%

# QV1 - (Fig. 5 upper half; in Elections paper)

QV1_AllRules <- filter(QV_AllRules, QV1_MinorityWon == TRUE) %>% 
  mutate(QV1_WelfareGain = ifelse(AggregateWelfare_QV1 > AggregateWelfare_noQV, TRUE, FALSE)) %>%
  group_by(Rule, QV1_WelfareGain) %>% summarize(Count = n()) %>% 
  mutate(Freq = round(Count/sum(Count), digits = 3)*100) %>% filter(QV1_WelfareGain == TRUE)

welfarePlot <- ggplot(data = QV1_AllRules, aes(x=Rule, y = Freq)) + geom_point() + 
  geom_text(aes(label=Freq), hjust = .5, vjust = 2) + ggtitle("QV1 (No Recal) - Freq. of Welfare Increases") + 
  scale_y_continuous(limits = c(0,100), breaks = seq(0, 100, by = 10)) +  
  xlab("Bonus Vote - Rules") + ylab("Freq. of Welfare Increases") + 
  theme(axis.text.x=element_text(size=12, vjust = -.05), legend.text=element_text(size=16), title=element_text(size=12))
print(welfarePlot)
ggsave(welfarePlot, file="./SV & QV Simulations (No Recalibration)/QV - Welfare & Minority Victories/QV1 (No Recal) - Freq. of Welfare Increases.png", height = 5.5, width = 7.5)

#QV2 - (Fig. 5 upper half; in Elections paper)

QV2_AllRules <- filter(QV_AllRules, QV2_MinorityWon == TRUE) %>% 
  mutate(QV2_WelfareGain = ifelse(AggregateWelfare_QV2 > AggregateWelfare_noQV, TRUE, FALSE)) %>%
  group_by(Rule, QV2_WelfareGain) %>% summarize(Count = n()) %>% 
  mutate(Freq = round(Count/sum(Count), digits = 3)*100) %>% filter(QV2_WelfareGain == TRUE)

welfarePlot <- ggplot(data = QV2_AllRules, aes(x=Rule, y = Freq)) + geom_point() + 
  geom_text(aes(label=Freq), hjust = .5, vjust = 2) + ggtitle("QV2 (No Recal) - Freq. of Welfare Increases") + 
  scale_y_continuous(limits = c(0,100), breaks = seq(0, 100, by = 10)) +  
  xlab("Bonus Vote - Rules") + ylab("Freq. of Welfare Increases") + 
  theme(axis.text.x=element_text(size=12, vjust = -.05), legend.text=element_text(size=16), title=element_text(size=12))
print(welfarePlot)
ggsave(welfarePlot, file="./SV & QV Simulations (No Recalibration)/QV - Welfare & Minority Victories/QV2 (No Recal) - Freq. of Welfare Increases.png", height = 5.5, width = 7.5)


# table of results (combined)

FreqWelfareIncrease <- bind_rows(mutate(QV1_AllRules, Vote = "QV1"), mutate(QV2_AllRules, Vote = "QV2")) %>% select(-matches("Welfare"))
sink("./SV & QV Simulations (No Recalibration)/QV - Welfare & Minority Victories/Table - Freq. of Welfare Gains (per Rule).txt")
stargazer(as.data.frame(FreqWelfareIncrease), type = "text", summary = FALSE, rownames = FALSE, digits = 1)
sink()


###################################################
####### Mean % Welfare Change Plot (byRule) #######
##################################################%

# QV1 -  (Fig. 5 bottom half; in Elections paper)
QV1_AllRules <- filter(QV_AllRules, QV1_MinorityWon == TRUE) %>% 
  mutate(QV1_PercentWelfareChange = (AggregateWelfare_QV1 - AggregateWelfare_noQV)/AggregateWelfare_noQV) %>%
  group_by(Rule) %>% summarize(QV1_MeanPercentWelfareChange = round(mean(QV1_PercentWelfareChange), digits = 3)*100) %>%
  rename(MeanPercentWelfareChange = QV1_MeanPercentWelfareChange)

welfarePlot <- ggplot(data = QV1_AllRules, aes(x=Rule, y = MeanPercentWelfareChange)) + 
  geom_bar(stat="identity") + scale_y_continuous(limits = c(-5,15), breaks = seq(-5, 15, by = 5)) +  
  xlab("Bonus Vote - Rules") + ylab("Mean Percent Welfare Change") + 
  theme(axis.text.x=element_text(size=12, vjust = -.05), legend.text=element_text(size=16), title=element_text(size=12)) +
  ggtitle("QV1 (No Recal) - Bootstrapped Mean Percent Welfare Change")
print(welfarePlot)
ggsave(welfarePlot, file="./SV & QV Simulations (No Recalibration)/QV - Welfare & Minority Victories/QV1 (No Recal) - Mean Percent Welfare Change.png", height = 5.5, width = 7.5)


# QV2 -  (Fig. 5 bottom half; in Elections paper)
QV2_AllRules <- filter(QV_AllRules, QV2_MinorityWon == TRUE) %>% 
  mutate(QV2_PercentWelfareChange = (AggregateWelfare_QV2 - AggregateWelfare_noQV)/AggregateWelfare_noQV) %>%
  group_by(Rule) %>% summarize(QV2_MeanPercentWelfareChange = round(mean(QV2_PercentWelfareChange), digits = 3)*100) %>%
  rename(MeanPercentWelfareChange = QV2_MeanPercentWelfareChange)

welfarePlot <- ggplot(data = QV2_AllRules, aes(x=Rule, y = MeanPercentWelfareChange)) + 
  geom_bar(stat="identity") + scale_y_continuous(limits = c(-5,15), breaks = seq(-5, 15, by = 5)) +  
  xlab("Bonus Vote - Rules") + ylab("Mean Percent Welfare Change") + 
  theme(axis.text.x=element_text(size=12, vjust = -.05), legend.text=element_text(size=16), title=element_text(size=12)) +
  ggtitle("QV2 (No Recal) - Bootstrapped Mean Percent Welfare Change")
print(welfarePlot)
ggsave(welfarePlot, file="./SV & QV Simulations (No Recalibration)/QV - Welfare & Minority Victories/QV2 (No Recal) - Mean Percent Welfare Change.png", height = 5.5, width = 7.5)


# table of results (combined)

MeanWelfareChange <- bind_rows(mutate(QV1_AllRules, Vote = "QV1"), mutate(QV2_AllRules, Vote = "QV2")) 
sink("./SV & QV Simulations (No Recalibration)/QV - Welfare & Minority Victories/Table - Mean Percent Welfare Change (per Rule).txt")
stargazer(as.data.frame(MeanWelfareChange), type = "text", summary = FALSE, rownames = FALSE, digits = 2)
sink()



###################################################
####### Hist of % Welfare Change (per Rule) #######
##################################################%

WelfareChangeHist <- function(sampleData, SetRule, SampleTitle)
{
  #--- Set Graphing Parameters
  
  bins_width5 <- c("[-20,-15)","[-15,-10)","[-10,-5)", "[-5,0)", "[0,5)", "[5,10)", "[10,15)", "[15,20)", "[20,25)", "[25,30)")
  bins5 <- data.frame(Bins = bins_width5, BinOrder = seq(1,length(bins_width5),1), stringsAsFactors = F)
  
  #--- Determines the current vote
  if(SampleTitle== "QV1"){
    tempRule <- filter(sampleData, QV1_MinorityWon == TRUE, Rule == SetRule) %>%
      mutate(QV1_PercentWelfareChange = (AggregateWelfare_QV1 - AggregateWelfare_noQV)/AggregateWelfare_noQV,
             QV1_PercentWelfareChange = round(QV1_PercentWelfareChange, digits = 3)*100,
             Bins = cut(QV1_PercentWelfareChange, right = F, breaks = seq(-20, 30, 5)))
  }
  
  if(SampleTitle == "QV2"){
    tempRule <- filter(sampleData, QV2_MinorityWon == TRUE, Rule == SetRule) %>%
      mutate(PercentWelfareChange = (AggregateWelfare_QV2 - AggregateWelfare_noQV)/AggregateWelfare_noQV,
             PercentWelfareChange = round(PercentWelfareChange, digits = 3)*100,
             Bins = cut(PercentWelfareChange, right = F, breaks = seq(-20, 30, 5))) 
  }
  
  tempRule <- tempRule %>% group_by(Bins) %>% summarize(Count = n()) %>% mutate(Freq = round(Count/sum(Count), digits = 3)*100) 
  tempRule <- full_join(bins5, tempRule, by = "Bins")
  tempRule[is.na(tempRule)] <- 0
  tempRule$Bins <- factor(tempRule$Bins, levels = tempRule$Bins[order(tempRule$BinOrder)])
  
  welfarePlot <- ggplot(tempRule, aes(x = Bins, y = Freq)) +
    geom_bar(stat="identity") + xlab("Percent Welfare Change") + ylab("Relative Frequency") +
    scale_y_continuous(limits = c(0,70), breaks = seq(0, 70, by = 10)) +
    theme(axis.text.x=element_text(size=12, vjust = .5, angle = 300), legend.text=element_text(size=16), title=element_text(size=12)) +
    ggtitle(paste("Rule ",SetRule, " (", SampleTitle, ") ", " - Histogram of % Welfare Change (No Recal)", sep = ""))
  print(welfarePlot)
  SetRule <- gsub("/",".", SetRule)
  # ggsave(welfarePlot, file=paste("./SV & QV Simulations (No Recalibration)/QV - Welfare & Minority Victories/Rule ", SetRule, " (", SampleTitle,") - Histogram of Percent Welfare Change (No Recal).png",sep=""), height = 4.5, width = 6.5)
  ggsave(welfarePlot, file=paste("./SV & QV Simulations (No Recalibration)/QV - Welfare & Minority Victories/Rule ", SetRule, " (", SampleTitle,") - Histogram of Percent Welfare Change.png",sep=""), height = 4.5, width = 6.5)
}

# QV1 - (Fig. 6 in Elections paper)

WelfareChangeHist(QV_AllRules, "A", "QV1")
WelfareChangeHist(QV_AllRules, "B", "QV1")
WelfareChangeHist(QV_AllRules, "B - CNV", "QV1")
WelfareChangeHist(QV_AllRules, "C", "QV1")
WelfareChangeHist(QV_AllRules, "D - Unif", "QV1")
WelfareChangeHist(QV_AllRules, "D - Emp", "QV1")
WelfareChangeHist(QV_AllRules, "D - 50/50", "QV1")
WelfareChangeHist(QV_AllRules, "E - Unif", "QV1")
WelfareChangeHist(QV_AllRules, "E - Emp", "QV1")

# QV2 - (Fig. 6 in Elections paper)

WelfareChangeHist(QV_AllRules, "A", "QV2")
WelfareChangeHist(QV_AllRules, "B", "QV2")
WelfareChangeHist(QV_AllRules, "B - CNV", "QV2")
WelfareChangeHist(QV_AllRules, "C", "QV2")
WelfareChangeHist(QV_AllRules, "D - Unif", "QV2")
WelfareChangeHist(QV_AllRules, "D - Emp", "QV2")
WelfareChangeHist(QV_AllRules, "D - 50/50", "QV2")
WelfareChangeHist(QV_AllRules, "E - Unif", "QV2")
WelfareChangeHist(QV_AllRules, "E - Emp", "QV2")





###########################################
##### Aggregate Preference Histograms #####
###########################################


ProposalHist <- function(sampleData, Proposal, RawPrefSummary)
{
  #--- Set Graphing Parameters
  
  plotMax <- 70
  
  bins_width10 <- c("[-100,-90)", "[-90,-80)", "[-80,-70)", "[-70,-60)", "[-60,-50)", "[-50,-40)", 
                    "[-40,-30)", "[-30,-20)", "[-20,-10)", "[-10,0)", "0", "(0,10]", "(10,20]", 
                    "(20,30]", "(30,40]", "(40,50]", "(50,60]", "(60,70]", "(70,80]", "(80,90]", "(90,100]")
  
  
  # Calculate the average total points for each of the proposals and store it on the dataframe
  PrefSummary <- summarise_all(RawPrefSummary, funs(mean))
  
  # Calculate the average number of subjects and store it on a temp dataframe
  temp <- summarise_all(sampleData, funs(mean))
  temp <- unname(unlist(temp[1,]))
  
  # Create dataframe with summarized aggregate data
  currentProp <- data.frame(Bins = factor(bins_width10), NumSubjects = temp, BinOrder = seq(1,21,1), 
                            BinColor = c(rep(0,10),1,rep(0,10)), stringsAsFactors = F)
  currentProp$Bins <- factor(currentProp$Bins, levels = currentProp$Bins[order(currentProp$BinOrder)])
  
  #plot
  intensityPlot <- ggplot(data = currentProp, aes(x=Bins, y = NumSubjects, fill=BinColor)) + 
    geom_bar(stat="identity") + scale_y_continuous(limits = c(0,plotMax), breaks = seq(0, plotMax, by = 10)) +  
    xlab("Number of Points Assigned to Proposal") + ylab("Number of Subjects") + 
    theme(axis.text.x=element_text(size=12, vjust = -.05, angle = 90), legend.text=element_text(size=16), title=element_text(size=14)) +
    ggtitle(paste("QV (No Recal) - ",Proposal, sep = "")) + guides(fill = F)
  
  l1 <- paste("Mean Total Points InFavor: ",PrefSummary[[paste(Proposal,"_PointsInFavor", sep = "")]], sep = "")
  l2 <- paste("Mean Total Points Against: ",PrefSummary[[paste(Proposal,"_PointsAgainst", sep = "")]], sep = "")
  intensityPlot <- intensityPlot + annotate("text", x = 18, y = c(70,66), label = c(l1,l2))
  
  print(intensityPlot)
  ggsave(intensityPlot, file=paste("SV & QV Simulations (No Recalibration)/QV - Preference Histograms/QV (No Recal) - ",Proposal,".png", sep = ""), height = 4.5, width = 9.5)
  
}

ProposalHist(sampleData = QV_ImmBSPref, Proposal = "Immigration", RawPrefSummary = QV_PrefSummary_raw)
ProposalHist(sampleData = QV_BondBSPref, Proposal = "PublicBonds", RawPrefSummary = QV_PrefSummary_raw)
ProposalHist(sampleData = QV_EducBSPref, Proposal = "BilingualEduc", RawPrefSummary = QV_PrefSummary_raw)
ProposalHist(sampleData = QV_TeachBSPref, Proposal = "TeacherTenure", RawPrefSummary = QV_PrefSummary_raw)


